#define _USE_MATH_DEFINES
#include <cmath>
#ifndef M_PI
#define M_PI 3.14159265358979323846
#endif

#include "data_structures.hpp"
#include <cmath>
#include <iostream>
#include <iomanip>
#include <numeric>
#include <limits>
#include <algorithm>

#ifdef _OPENMP
#include <omp.h>
#endif

namespace EddyDetection {
namespace Utils {

VelocityField createDoubleGyreField(int grid_x, int grid_y, double t) {
    // Double eddy parameters
    const double A = 0.1;
    const double epsilon = 0.25;
    const double omega = 2.0 * M_PI / 10.0;
    
    // Create coordinate grid
    Eigen::VectorXd lat = Eigen::VectorXd::LinSpaced(grid_x, 0.0, 2.0);
    Eigen::VectorXd lon = Eigen::VectorXd::LinSpaced(grid_y, 0.0, 1.0);
    
    // Initialize velocity field
    Eigen::MatrixXd u(grid_x, grid_y);
    Eigen::MatrixXd v(grid_x, grid_y);
    
    // Calculate time related parameters
    double a = epsilon * std::sin(omega * t);
    double b = 1.0 - 2.0 * epsilon * std::sin(omega * t);
    
    for (int i = 0; i < grid_x; ++i) {
        for (int j = 0; j < grid_y; ++j) {
            double x = lat[i];
            double y = lon[j];
            double f = a * x * x + b * x;
            double df_dx = 2.0 * a * x + b;
            
            // Calculate the velocity component
            u(i, j) = -M_PI * A * std::sin(M_PI * f) * std::cos(M_PI * y);
            v(i, j) = M_PI * A * std::cos(M_PI * f) * std::sin(M_PI * y) * df_dx;
        }
    }
    
    VelocityField field(u, v, lat, lon);
    field.time = Eigen::VectorXd::Constant(1, t); // Single time point
    
    return field;
}

bool validateVelocityField(const VelocityField& field) {

    if (field.u_layers.empty() || field.v_layers.empty()) {
        std::cerr << "Error: Speed field data is empty" << std::endl;
        return false;
    }
    
    const auto& u = field.u();
    const auto& v = field.v();
    
    if (u.rows() != field.lat.size()) {
        std::cerr << "Error: The number of rows in field u does not match the number of latitude coordinates" << std::endl;
        return false;
    }
    
    if (u.cols() != field.lon.size()) {
        std::cerr << "Error: The number of columns in field u does not match the number of longitude coordinates" << std::endl;
        return false;
    }
    
    if (v.rows() != u.rows() || v.cols() != u.cols()) {
        std::cerr << "Error: u and v velocity field dimensions do not match" << std::endl;
        return false;
    }
    
    for (int i = 1; i < field.lat.size(); ++i) {
        if (field.lat[i] <= field.lat[i-1]) {
            std::cerr << "Warning: Latitude coordinates are not strictly monotonically increasing" << std::endl;
            break;
        }
    }
    
    for (int i = 1; i < field.lon.size(); ++i) {
        if (field.lon[i] <= field.lon[i-1]) {
            std::cerr << "Warning: Longitude coordinates are not strictly monotonically increasing" << std::endl;
            break;
        }
    }
    

    int nan_count = 0;
    int inf_count = 0;
    
    for (int i = 0; i < u.rows(); ++i) {
        for (int j = 0; j < u.cols(); ++j) {
            if (std::isnan(u(i, j)) || std::isnan(v(i, j))) {
                nan_count++;
            }
            if (std::isinf(u(i, j)) || std::isinf(v(i, j))) {
                inf_count++;
            }
        }
    }
    
    if (nan_count > 0) {
        std::cout << "Information: The velocity field includes " << nan_count << " NaN value" << std::endl;
    }
    
    if (inf_count > 0) {
        std::cerr << "Warning: The speed field contains " << inf_count << " Infinite value" << std::endl;
    }
    
    return true;
}

bool validateTemperatureField(const TemperatureField& field) {

    if (field.data_layers.empty()) {
        std::cerr << "Error: Temperature field data is empty" << std::endl;
        return false;
    }
    
    const auto& data = field.data();
    
    if (data.rows() != field.lat.size()) {
        std::cerr << "Error: The number of rows in the temperature field does not match the number of latitude coordinates" << std::endl;
        return false;
    }
    
    if (data.cols() != field.lon.size()) {
        std::cerr << "Error: The number of columns in the temperature field does not match the number of longitude coordinates" << std::endl;
        return false;
    }

    for (int i = 1; i < field.lat.size(); ++i) {
        if (field.lat[i] <= field.lat[i - 1]) {
            std::cerr << "Warning: Latitude coordinates are not strictly monotonically increasing" << std::endl;
            break;
        }
    }

    for (int i = 1; i < field.lon.size(); ++i) {
        if (field.lon[i] <= field.lon[i-1]) {
            std::cerr << "Warning: Longitude coordinates are not strictly monotonically increasing" << std::endl;
            break;
        }
    }
     

    double min_temp = data.minCoeff();
    double max_temp = data.maxCoeff();
    
    if (min_temp < -10.0 || max_temp > 50.0) {
        std::cerr << "Warning: Temperature range may not be reasonable (" << min_temp << "°C to " << max_temp << "°C)" << std::endl;
    }
    

    int nan_count = 0;
    int inf_count = 0;
    
    for (int i = 0; i < data.rows(); ++i) {
        for (int j = 0; j < data.cols(); ++j) {
            if (std::isnan(data(i, j))) {
                nan_count++;
            }
            if (std::isinf(data(i, j))) {
                inf_count++;
            }
        }
    }
    
    if (nan_count > 0) {
        std::cout << "Information: The temperature field includes " << nan_count << " NaN value" << std::endl;
    }
    
    if (inf_count > 0) {
        std::cerr << "Warning: The temperature field contains " << inf_count << " Infinite value" << std::endl;
    }
    
    return true;
}

double degToKm(double deg_diff, double mean_lat_rad) {
    const double R_EARTH = 6371.0;  // Earth radius (kilometers)
    const double deg_to_rad = M_PI / 180.0;
    
    // Convert longitude difference to kilometers
    double km_per_deg_lon = R_EARTH * deg_to_rad * std::cos(mean_lat_rad);
    
    return deg_diff * km_per_deg_lon;
}

Eigen::Vector2d pixelToCoordinate(const Eigen::Vector2i& pixel,
                                  const Eigen::VectorXd& lat,
                                  const Eigen::VectorXd& lon) {
    if (pixel[0] < 0 || pixel[0] >= lat.size() || 
        pixel[1] < 0 || pixel[1] >= lon.size()) {
        return Eigen::Vector2d(0.0, 0.0);
    }
    
    return Eigen::Vector2d(lat[pixel[0]], lon[pixel[1]]);
}

EddyDetectionConfig suggestOptimalParameters(const VelocityField& velocity_field) {
    EddyDetectionConfig config;
    
    const auto& lat = velocity_field.lat;
    const auto& lon = velocity_field.lon;
    const auto& u = velocity_field.u();
    const auto& v = velocity_field.v();
    
    std::cout << "\n=== Data feature analysis and parameter optimization suggestions ===" << std::endl;
    
    if (u.size() == 0 || v.size() == 0) {
        std::cerr << "Error: Speed field data is empty" << std::endl;
        return config;
    }
    
    double mean_lat_rad = lat.mean() * M_PI / 180.0;
    double dy_km = (lat.size() > 1 ? (lat[1] - lat[0]) : 1.0) * 111.32;
    double dx_km = (lon.size() > 1 ? (lon[1] - lon[0]) : 1.0) * 111.32 * std::cos(mean_lat_rad);
    double resolution_km = (dx_km + dy_km) / 2.0;
    
    std::cout << "grid resolution: ~" << std::fixed << std::setprecision(1) << resolution_km << " km" << std::endl;
    
    double suggested_min_area = std::pow(resolution_km * 5, 2);
    config.min_area_km2 = suggested_min_area;
    std::cout << "Suggested minimum Eddy area: " << std::fixed << std::setprecision(0) << suggested_min_area << " km²" << std::endl;
    
    // Calculate the average flow velocity
    std::vector<double> speed_values;
    for (int i = 0; i < u.rows(); ++i) {
        for (int j = 0; j < u.cols(); ++j) {
            if (std::isfinite(u(i, j)) && std::isfinite(v(i, j))) {
                double speed = std::sqrt(u(i, j) * u(i, j) + v(i, j) * v(i, j));
                speed_values.push_back(speed);
            }
        }
    }
    
    if (!speed_values.empty()) {
        double mean_speed = std::accumulate(speed_values.begin(), speed_values.end(), 0.0) / speed_values.size();
        std::cout << "average velocity: " << std::fixed << std::setprecision(3) << mean_speed << " m/s" << std::endl;
        
        double suggested_distance = std::sqrt(suggested_min_area);
        double suggested_time_hours = (suggested_distance * 1000) / mean_speed / 3600;
        config.integration_time_hours = std::max(24.0, std::min(120.0, suggested_time_hours));
        std::cout << "Suggested points time: " << std::fixed << std::setprecision(0) << config.integration_time_hours 
                  << " hour (" << config.integration_time_hours/24.0 << " day)" << std::endl;
    }
    
    if (resolution_km < 5.0) {
        config.ftle_threshold_percentile = 97.0;
        config.vorticity_threshold_percentile = 90.0;
        std::cout << "High resolution data detection using strict thresholds" << std::endl;
    } else if (resolution_km > 20.0) {
        config.ftle_threshold_percentile = 90.0;
        config.vorticity_threshold_percentile = 80.0;
        std::cout << "Low resolution data detection, using a looser threshold" << std::endl;
    } else {
        config.ftle_threshold_percentile = 95.0;
        config.vorticity_threshold_percentile = 87.5;
        std::cout << "Medium resolution data detection, using optimized default threshold values" << std::endl;
    }
    
    std::cout << "Final suggested parameters:" << std::endl;
    std::cout << "  - FTLE Threshold: " << config.ftle_threshold_percentile << "%" << std::endl;
    std::cout << "  - Vorticity Threshold: " << config.vorticity_threshold_percentile << "%" << std::endl;
    std::cout << "  - Minimum Area: " << config.min_area_km2 << " km²" << std::endl;
    std::cout << "  - Integration Time: " << config.integration_time_hours << " hour" << std::endl;
    
    return config;
}

void printDataCharacteristics(const VelocityField& velocity_field) {
    const auto& lat = velocity_field.lat;
    const auto& lon = velocity_field.lon;
    const auto& u = velocity_field.u();
    const auto& v = velocity_field.v();
    
    if (u.size() == 0 || v.size() == 0) {
        std::cerr << "Error: Speed field data is empty" << std::endl;
        return;
    }
    
    std::cout << "\n=== Data consistency check ===" << std::endl;
    std::cout << "Latitude range: " << std::fixed << std::setprecision(2)
        << lat.minCoeff() << "° to " << lat.maxCoeff() << "°" << std::endl;
    std::cout << "Longitude range: " << std::fixed << std::setprecision(2) 
              << lon.minCoeff() << "° to " << lon.maxCoeff() << "°" << std::endl;

    if (lat.size() > 1) {
        double mean_dlat = (lat[lat.size() - 1] - lat[0]) / (lat.size() - 1);
        std::cout << "Latitude step size: " << std::fixed << std::setprecision(4) << mean_dlat << "°" << std::endl;
    }

    if (lon.size() > 1) {
        double mean_dlon = (lon[lon.size()-1] - lon[0]) / (lon.size() - 1);
        std::cout << "Longitude step size: " << std::fixed << std::setprecision(4) << mean_dlon << "°" << std::endl;
    }
    
    std::cout << "U velocity field shape: " << u.rows() << " × " << u.cols() << std::endl;
    std::cout << "V velocity field shape: " << v.rows() << " × " << v.cols() << std::endl;
    
    int total_points = u.rows() * u.cols();
    int valid_u = 0, valid_v = 0;
    double u_min = std::numeric_limits<double>::max();
    double u_max = std::numeric_limits<double>::lowest();
    double v_min = std::numeric_limits<double>::max();
    double v_max = std::numeric_limits<double>::lowest();
    
    for (int i = 0; i < u.rows(); ++i) {
        for (int j = 0; j < u.cols(); ++j) {
            if (std::isfinite(u(i, j))) {
                valid_u++;
                u_min = std::min(u_min, u(i, j));
                u_max = std::max(u_max, u(i, j));
            }
            if (std::isfinite(v(i, j))) {
                valid_v++;
                v_min = std::min(v_min, v(i, j));
                v_max = std::max(v_max, v(i, j));
            }
        }
    }
    
    std::cout << "Effective value ratio of U speed: " << std::fixed << std::setprecision(1) 
              << (100.0 * valid_u / total_points) << "%" << std::endl;
    std::cout << "Effective value ratio of V speed: " << std::fixed << std::setprecision(1) 
              << (100.0 * valid_v / total_points) << "%" << std::endl;
    
    if (valid_u > 0) {
        std::cout << "U speed range: " << std::fixed << std::setprecision(3) 
                  << u_min << " to " << u_max << " m/s" << std::endl;
    }
    
    if (valid_v > 0) {
        std::cout << "V speed range: " << std::fixed << std::setprecision(3) 
                  << v_min << " to " << v_max << " m/s" << std::endl;
    }
    
    if (valid_u > 0 && valid_v > 0) {
        std::vector<double> speed_values;
        for (int i = 0; i < u.rows(); ++i) {
            for (int j = 0; j < u.cols(); ++j) {
                if (std::isfinite(u(i, j)) && std::isfinite(v(i, j))) {
                    double speed = std::sqrt(u(i, j) * u(i, j) + v(i, j) * v(i, j));
                    speed_values.push_back(speed);
                }
            }
        }
        
        if (!speed_values.empty()) {
            double mean_speed = std::accumulate(speed_values.begin(), speed_values.end(), 0.0) / speed_values.size();
            std::cout << "Speed magnitude: " << std::fixed << std::setprecision(3) << mean_speed << " m/s" << std::endl;
        }
    }
}

} // namespace Utils

// Main interface class implementation
EddyDetector::EddyDetector(const EddyDetectionConfig& config) : config_(config) {
    std::cout << "Eddy detector initialization completed" << std::endl;
    std::cout << "Configuration Parameter: " << std::endl;
    std::cout << "  - FTLE threshold percentile: " << config_.ftle_threshold_percentile << "%" << std::endl;
    std::cout << "  - Vorticity threshold percentile: " << config_.vorticity_threshold_percentile << "%" << std::endl;
    std::cout << "  - Minimum Area: " << config_.min_area_km2 << " km²" << std::endl;
    std::cout << "  - integration time: " << config_.integration_time_hours << " hour" << std::endl;
    std::cout << "  - Use OpenMP: " << (config_.use_openmp ? "Yes" : "No") << std::endl;
}

std::string EddyDetector::detectEddiesLagrangian(const VelocityField& velocity_field,
                                                 const TemperatureField& temp_field,
                                                 double depth_m,
                                                 int time_index) {
    std::cout << "\n=== Start Lagrange Eddy detection ===" << std::endl;
    
    if (!Utils::validateVelocityField(velocity_field)) {
        std::cerr << "Speed field verification failed" << std::endl;
        return "";
    }
    
    if (!Utils::validateTemperatureField(temp_field)) {
        std::cerr << "Temperature field verification failed" << std::endl;
        return "";
    }
    
    Utils::printDataCharacteristics(velocity_field);
    EddyDetectionConfig optimal_config = Utils::suggestOptimalParameters(velocity_field);
    
    std::cout << "\n Using optimized parameters for detection..." << std::endl;
    
    auto eddies = LagrangianCore::detectEddiesLagrangian(velocity_field, temp_field, optimal_config, time_index);
    
    std::string geojson = GeoJSONExporter::exportToGeoJSON(eddies, "Lagrangian", depth_m);
    
    std::cout << "Using optimized parameters for detection, Lagrange method detection was completed, and a total of findings were found " 
        << eddies.size() << " eddy" << std::endl;
    return geojson;
}

std::string EddyDetector::detectEddiesEulerian(const VelocityField& velocity_field,
                                               const TemperatureField& temp_field,
                                               double depth_m,
                                               int time_index) {
    std::cout << "\n=== Start Euler Eddy detection ===" << std::endl;
    
    if (!Utils::validateVelocityField(velocity_field)) {
        std::cerr << "Speed field verification failed" << std::endl;
        return "";
    }
    
    if (!Utils::validateTemperatureField(temp_field)) {
        std::cerr << "Temperature field verification failed" << std::endl;
        return "";
    }
    
    Utils::printDataCharacteristics(velocity_field);
    EddyDetectionConfig optimal_config = Utils::suggestOptimalParameters(velocity_field);
    
    std::cout << "\n Using optimized parameters for detection..." << std::endl;
    
    auto eddies = EulerianCore::detectEddiesEulerian(velocity_field, temp_field, optimal_config, time_index);
    
    std::string geojson = GeoJSONExporter::exportToGeoJSON(eddies, "Eulerian", depth_m);
    
    std::cout << "Euler method detection completed, found a total of " << eddies.size() << " Eddy" << std::endl;
    return geojson;
}

std::pair<std::string, std::string> EddyDetector::detectEddiesBothMethods(
    const VelocityField& velocity_field,
    const TemperatureField& temp_field,
    double depth_m,
    int time_index) {
    
    std::cout << "\n=== Start dual method Eddy detection ===" << std::endl;
    
    std::string lagrangian_geojson = detectEddiesLagrangian(velocity_field, temp_field, depth_m, time_index);
    std::string eulerian_geojson = detectEddiesEulerian(velocity_field, temp_field, depth_m, time_index);
    
    std::cout << "Dual method detection completed" << std::endl;
    
    return std::make_pair(lagrangian_geojson, eulerian_geojson);
}




// =============================== START With NO Lingking 3D =====================================
std::vector<std::pair<double, std::vector<EddyContourData>>> Eddy::star2D(const VelocityData& velocity_data,
    const TemperatureData& temp_data,
    const EddyDetectionConfig& config,
    const std::string& method) {
    VelocityField velocity_field;
    TemperatureField temp_field;

    DataTransform transform;
    transform.data_to_field(velocity_data, temp_data, velocity_field, temp_field);

    EddyDetector detection(config);
    std::vector<std::vector<EddyContour>> all_layers_eddies = detection.detectLayers2DEddies(velocity_field, temp_field, method);

    std::vector<std::vector<EddyContourData>> all_layers_eddiesdata;
    transform.eddy2D_to_2DData(all_layers_eddiesdata, all_layers_eddies);

    std::vector<std::pair<double, std::vector<EddyContourData>>> result;
    for (size_t i = 0; i < all_layers_eddiesdata.size(); ++i) {
        double depth = velocity_field.depth[i];
        result.emplace_back(depth, all_layers_eddiesdata[i]);
    }
    return result;
}

std::vector<std::vector<EddyContour>> EddyDetector::detectLayers2DEddies(const VelocityField& velocity_field,
    const TemperatureField& temp_field,
    const std::string& method) {
    std::cout << "\n=== Start 2D Eddy detection (" << method << " method) ===" << std::endl;

    size_t num_layers = velocity_field.depth.size();
    if (num_layers == 0) {
        std::cerr << "Error: There are no deep layers in the data." << std::endl;
        return std::vector<std::vector<EddyContour>>();
    }
    std::cout << "We need to analyze a total of  " << num_layers << " deep layer." << std::endl;
    size_t num_u_layers = velocity_field.u_layers.size();
    std::cout << "the size of u is " << num_u_layers << std::endl;

    std::vector<std::vector<EddyContour>> all_layers_eddies(num_layers);

    // 1. Parallel layer by layer 2D Eddy detection
    #pragma omp parallel for if(config_.use_openmp)
    for (int i = 0; i < num_layers; ++i) {
        const auto& u = velocity_field.u_layers[i];
        const auto& v = velocity_field.v_layers[i];
        const auto& temp = temp_field.data_layers[i];
        
        auto eddies_2d = detect2DEddiesForLayer(u, v, temp, velocity_field.lat, velocity_field.lon, method);
        
        for(auto& eddy : eddies_2d) {
            eddy.depth = velocity_field.depth[i];
            eddy.layer_index = i;
        }

        all_layers_eddies[i] = eddies_2d;
    
        #pragma omp critical
        {
            if (!eddies_2d.empty()) {
                std::cout << "  -> depth " << velocity_field.depth[i] << "m (layer " << i << "): Detected " << eddies_2d.size() << " Eddy." << std::endl;
            }
        }
    }

    //// 1. Sequential layer by layer 2D Eddy detection
    //for (int i = 0; i < num_layers; ++i) {
    //    const auto& u = velocity_field.u_layers[i];
    //    const auto& v = velocity_field.v_layers[i];
    //    const auto& temp = temp_field.data_layers[i];
    //
    //    auto eddies_2d = detect2DEddiesForLayer(u, v, temp, velocity_field.lat, velocity_field.lon, method);
    //
    //    for (auto& eddy : eddies_2d) {
    //        eddy.depth = velocity_field.depth[i];
    //        eddy.layer_index = i;
    //    }
    //
    //    all_layers_eddies[i] = eddies_2d;
    //
    //    if (!eddies_2d.empty()) {
    //        std::cout << "  -> depth " << velocity_field.depth[i] << "m (layer " << i << "): Detected " << eddies_2d.size() << " Eddy." << std::endl;
    //    }
    //}
    return all_layers_eddies;
}
// =============================== END With NO Lingking 3D =====================================
// =============================== START With Lingking 3D =====================================
std::vector<Eddy3DData> Eddy::star(const VelocityData& velocity_data,
                                        const TemperatureData& temp_data,
                                    const EddyDetectionConfig& config,
                                    const std::string& method) {
    VelocityField velocity_field;
    TemperatureField temp_field;
    
    DataTransform transform;
    transform.data_to_field(velocity_data, temp_data, velocity_field, temp_field);
    
    EddyDetector detection(config);
    std::vector<Eddy3D> eddies_3d = detection.detect3DEddies(velocity_field, temp_field, method);
    
    std::vector<Eddy3DData> eddiesdata_3d;
    transform.eddies_to_eddiesdata(eddiesdata_3d, eddies_3d);
    
    return eddiesdata_3d;
}

std::vector<Eddy3D> EddyDetector::detect3DEddies(const VelocityField& velocity_field,
                                       const TemperatureField& temp_field,
                                       const std::string& method) {
    std::cout << "\n=== Start 3D Eddy detection (" << method << " method) ===" << std::endl;
    
    size_t num_layers = velocity_field.depth.size();
    if (num_layers == 0) {
        std::cerr << "Error: There are no deep layers in the data." << std::endl;
        return std::vector<Eddy3D>();
    }
    std::cout << "We need to analyze a total of  " << num_layers << " deep layer." << std::endl;
    size_t num_u_layers = velocity_field.u_layers.size();
    std::cout << "the size of u is " << num_u_layers << std::endl;
    
    std::vector<std::vector<EddyContour>> all_layers_eddies(num_layers);

    //// 1. Parallel layer by layer 2D Eddy detection
    //#pragma omp parallel for if(config_.use_openmp)
    //for (int i = 0; i < num_layers; ++i) {
    //    const auto& u = velocity_field.u_layers[i];
    //    const auto& v = velocity_field.v_layers[i];
    //    const auto& temp = temp_field.data_layers[i];
    //    
    //    auto eddies_2d = detect2DEddiesForLayer(u, v, temp, velocity_field.lon, velocity_field.lat, method);
    //    
    //    for(auto& eddy : eddies_2d) {
    //        eddy.depth = velocity_field.depth[i];
    //        eddy.layer_index = i;
    //    }
    //    
    //    all_layers_eddies[i] = eddies_2d;
    //
    //    #pragma omp critical
    //    {
    //        if (!eddies_2d.empty()) {
    //            std::cout << "  -> depth " << velocity_field.depth[i] << "m (layer " << i << "): Detected " << eddies_2d.size() << " Eddy." << std::endl;
    //        }
    //    }
    //}

    // 1. Sequential layer by layer 2D Eddy detection
    for (int i = 0; i < num_layers; ++i) {
        const auto& u = velocity_field.u_layers[i];
        const auto& v = velocity_field.v_layers[i];
        const auto& temp = temp_field.data_layers[i];

        auto eddies_2d = detect2DEddiesForLayer(u, v, temp, velocity_field.lat, velocity_field.lon, method);

        for (auto& eddy : eddies_2d) {
            eddy.depth = velocity_field.depth[i];
            eddy.layer_index = i;
        }

        all_layers_eddies[i] = eddies_2d;

        if (!eddies_2d.empty()) {
            std::cout << "  -> depth " << velocity_field.depth[i] << "m (layer " << i << "): Detected " << eddies_2d.size() << " Eddy." << std::endl;
        }
    }
    
    // 2. Linking 3D structures
    std::cout << "\nLinking 3D structures..." << std::endl;
    std::vector<Eddy3D> eddies_3d = ThreeDLinker::link3DEddies(all_layers_eddies, config_);
    std::cout << "  -> Constructed " << eddies_3d.size() << " three-dimensional Eddy structures" << std::endl;

    return eddies_3d;
}
// =============================== END With Lingking 3D =====================================

std::vector<EddyContour> EddyDetector::detect2DEddiesForLayer(
    const Eigen::MatrixXd& u, const Eigen::MatrixXd& v, 
    const Eigen::MatrixXd& temp, const Eigen::VectorXd& lat,
    const Eigen::VectorXd& lon, const std::string& method) {
        
    std::vector<EddyContour> eddies;
    if (method == "Eulerian") {
        // Create temporary Velocity Field and Temperature Field
        VelocityField temp_vel_field(u, v, lat, lon);
        std::vector<Eigen::MatrixXd> temp_layers;
        temp_layers.push_back(temp);
        TemperatureField temp_temp_field(temp_layers, lat, lon);
        
        // Using Euler's method for detection
        eddies = EulerianCore::detectEddiesEulerian(temp_vel_field, temp_temp_field, config_, 0);
    } else if (method == "Lagrangian") {
        // Create temporary Velocity Field and Temperature Field
        VelocityField temp_vel_field(u, v, lat, lon);
        std::vector<Eigen::MatrixXd> temp_layers;
        temp_layers.push_back(temp);
        TemperatureField temp_temp_field(temp_layers, lat, lon);
        
        // Detecting using Lagrange method
        eddies = LagrangianCore::detectEddiesLagrangian(temp_vel_field, temp_temp_field, config_, 0);
    }
    
    // Assign temporary IDs to the vortices in this layer
    for (size_t i = 0; i < eddies.size(); ++i) {
        eddies[i].id = "2D-" + std::to_string((int)eddies[i].type) + "-" + std::to_string(i);
    }
    
    return eddies;
}

} // namespace EddyDetection 